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Small localized black holes in a braneworld: 
Formulation and numerical method 
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No realistic black holes localized on a 3-brane in the Randall-Sundrum infinite braneworld have 
been found so far. The problem of finding a static black hole solution is reduced to a boundary value 
problem. We solve it by means of a numerical method, and show numerical examples of a localized 
black hole whose horizon radius is small compared to the bulk curvature scale. The sequence of 
small localized black holes exhibits a smooth transition from a five-dimensional Schwarzschild black 
hole, which is a solution in the limit of small horizon radius. The localized black hole tends to 
fiatten as its horizon radius increases. However, it becomes difficult to find black hole solutions as 
its horizon radius increases. 

PACS numbers; 04.50.+h, 04.70.Bw, 04.70.Dy 



I. INTRODUCTION 



. . . 

' Higher-dimensional black holes have been considered for a long time as purely theoretical applications motivated 
by higher-dimensional theories, such as string theory. However, recent developments in the scenario of large extra 
CO dimensions 1] have aroused new interest in such black holes. In the braneworld scenario, an interesting possibility 
H^' ] of black hole production at a collider was pointed out 0, 0] (see, e.g^ H and references therein) . There are other 
. types of braneworld models proposed by Randall and Sundrum (RS) [3,1a- In these models, the geometry warped in 
the direction of an extra dimension is used to explain the hierarchy between the TeV scale and the Planck scale, and 
to realize four-dimensional gravity effectively on the 3-brane. Also in the context of RS models, higher-dimensional 
' black holes may play an important roll. These so-called braneworld scenarios provide new and interesting situations 
CO ', to investigate higher-dimensional black holes. 

■ In the model of large extra dimensions, a physically meaningful sequence of black hole solutions will be obtained 
I as a slight modification of the higher-dimensional Kerr black hole (or more simply the Schwarzschild black hole 
O^' 0)S^3) if the horizon radius is sufficiently small compared to the extension of extra dimensions and the self-gravity 
^ ' due to periodic boundary is weak. As another sequence, there are the black string solutions. A black string is in 
(2J[) general unstable to linear perturbations with long wavelength in the direction along the string, which is called Gregory- 
Laflamme instability • (^ stability analysis of black strings is also found in Ref . [13.) Therefore, a black string 
is unstable when the horizon radius of the black string is sufficiently small compared to the extension of the extra 
dimension. For such a small black hole, the former sequence is expected to be stable. For RS models, since the 3-brane 
J-j I has tension, it is more difficult to find black hole solutions. Trivial black string solution is allowed also in these models, 
and it becomes unstable in the same way. There are many discussions about black holes in this model and some black 
hole solutions have been considered by several authors [ll[ll|ll[lllIl[ll[ll|2l[2l|2l[2l[M|2llME^ A 
strategy to construct a black hole solution is to assume an induced metric on the 3-brane as initial data, and extend 
it to the bulk analytically or numerically. This method generally results in a naked singularity in the bulk since there 
is no guarantee that the induced metric assumed as a boundary condition is compatible with a regular geometry. If 
we randomly specify the boundary metric on the 3-brane, almost all solutions develop a naked singularity. After all, 
no realistic black hole solutions which are stable and have no naked singularity have been found so far, and finding 
them is an interesting open question of nonlinear gravity in the braneworld. 

Successful recovery mechanism of 4D gravity on the 3-brane [2^ [s^, IS 112, IH, H H^l suggests the existence 

of black holes in the RS infinite braneworld f^. The shape of black holes is conjectured based on the Gregory-Laflamme 
instability (26l |. It was argued that an unstable black string will be pinched into many black holes. In addition to 
black holes in the bulk, there will be a black hole that is localized on the 3-brane. We can imagine this as a black hole 
bound to a domain wall if the 3-brane is realized by a domain wall |38|| . An exact solution representing a localized 
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black hole is known in the 4D braneworld model, which we call in this paper the Emparan-Horowitz-Myers (EHM) 
solution ^39]. However the corresponding solution in the original 5D braneworld model has not been discovered. While 
the localization of a black hole was motivated by the classical instability of a black string, the anticipated dynamics 
of pinching off the horizon is questioned by Horowitz and Maeda H(f\. Recently, Wiseman discovered nonuniform 
black strings by numerical calculation (4lj. The obtained solutions are likely to be unstable, but they are suggestive 
of a missing link between the localized black hole and the black string. On the other hand, based on an extensive 
use of the AdS/CFT correspondence, there is another discussion anticipating the absence of static localized black 
holes 42, 43]. 

The present paper explores the problem of black holes in the RS infinite braneworld. We consider a numerical 
construction of a black hole solution in this model. The method we use in this paper is based on a scheme developed 
by Wiseman |4^ . Our method does not require any assumptions for the induced metric on the 3-brane. We solve 
Einstein equations numerically under the boundary conditions determined by physical requirements. We will find 
small black hole solutions, whose horizon radius is smaller than the AdS curvature radius, although large black hole 
solutions have not been discovered. 

In the next section, we explain our method formulating the problem to be suitable for numerical calculation, and 
the boundary conditions are also discussed. The examples of numerical solutions are shown in Sec. IIIII Section [l VI is 
devoted to a summary. 



II. FORMULATION AS A BOUNDARY VALUE PROBLEM 



A. Conformal transformation 



To obtain black hole solutions in the RS infinite braneworld model, it is crucial to formulate the problem as a 
boundary value problem without assum ing any artificial boundary conditions. We follow and develop the numerical 
method developed in Ref. 13] (see also What we want to find is a static black hole solution that is localized on 

the brane. For simplicity, we further restrict our attention to the configuration with 5D axial symmetry (4D spherical 
symmetry). Under the static and axisymmetric assumption, the metric depends only on the radial coordinate r and 
the coordinate z in the direction of the extra dimension. Then, without loss of generality, the metric can be written 
as 

ds^ = ^ [-T^dt^ + e2«(dr2 + dz^) + r^e^^dfl^] , (1) 

where dfi^ — dO^ + sin^ Odcj)^ represents the line element on a unit 2-sphere. The cosmological constant in the bulk is 
related to the bulk curvature length £ by A = — 6/£^. If we set T = 1 and i? = C = 0, this metric becomes the AdS 
metric in the Poincare coordinates. Since we consider a localized black hole, polar coordinates 

r psinx, 
z — £ + pcosxj 

^ = X\ (2) 

are more convenient. Then, we have 

2 

dr^ + dz^ = dp^ + ^d^2_ (3) 

The angular coordinate ^ is useful to treat the coordinate singularity at x = numerically since the singularity 
becomes milder in ^ coordinate as we will see later [45l |. 

For numerical calculations, it is convenient if the boundaries are located on lines where one of the coordinates is 
constant. This can be achieved in general by using the residual gauge degrees of freedom. The metric form has 
the gauge degrees of freedom of conformal transformations in the two-dimensional space spanned by r and z. For 
our present purpose, it is convenient to use the conformal polar coordinates (C,x)j where C, = log p. Using these 
coordinates, we have dr^ + dz^ — p^idC,"^ + dx^)- For the conformal transformation C' = /(C; x) and x' = diCi X)^ f 
and g satisfy the Cauchy-Riemann relations. With a function 4'iCiX)j we have 

dc'+dx'^ = ^/j{c,x){de+dx'). (4) 



We can use these conformal degrees of freedom to deform the shape of the boundaries. From the Cauchy-Riemann 
relations, g (and also /) satisfies Laplace equation, and thus g is completely determined by specifying its boundary 
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conditions. Suppose we set g = on the axis and g = 7r/2 on the brane, and impose Neumann boundary conditions on 
the horizon boundary and the asymptotic boundary. The transformation / is determined by integrating the Cauchy- 
Riemann relations. The Neumann conditions for g on the horizon boundary and on the asymptotic boundary guarantee 
the constancy of / on these two boundaries. In this manner, we can generally find a conformal transformation 
that transforms each boundary to a constant coordinate line. Therefore, the location of the event horizon can be 
transformed to be 

Plhorizon = const Ph , (5) 

and the brane can be placed at z = £ |30j |. Here we note that / is uniquely determined up to adding a constant, and 
hence, the difference of the value of between two points cannot be changed arbitrarily. We therefore do not have a 
degree of freedom to change the ratio p/ ph for a given solution. 



B. Elliptic equations and constraint equations 



Let us consider the elliptic equations and constraint equations that we solve. From the five dimensional vacuum 
Einstein equations, C/^ :— — |Ag^ = 0, we obtain an elliptic equation for each metric component. Respective 

equations are given from t/*, {Ql — Gp — + "^Gg) and Gg as 
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where a comma means a partial derivative, and A = d^. + = dp + p ^dp + 2p ^(2^9| + d^) is the Laplace operator. 
The first terms in the equations for R and C have the factor that behaves as cx x~^ when x 0. Introduction of 
^ changes this severe behavior and it makes the terms more tractable in numerical calculation. Constraint equations 
are obtained, respectively, from t/| and {Gl — Gp + Gf + 2Gg) as 
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We perform numerical calculations for the elliptic equations using a relaxation method. The constraint equations 
are not explicitly solved but they are used as a check of accuracy. Before proceeding further, we discuss an important 
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property of the constraint equations |4J| . The constraint equations with an appropriate function multiphed satisfy 
Cauchy-Riemann relations. They are obtained from the nontrivial components of the Bianchi identities V^^/'^^ = 0. 

Assuming that equations for T, R, and C are satisfied, i.e., Gt = Qg = {Gp + S|) = 0, we obtain Cauchy-Riemann 
relations for t/| and {G^ ^ ^|) 

dcU~d^V = 0, 

d^U + d(^V = 0, (11) 
where U :== ■= PV^V^iG^ - e^|)/(4sin6l), and g := det{g^^). Then 

v= 0)'£!;!M!e,. ,12) 

Since U and V satisfy Cauchy-Riemann relations, each of them satisfies the Laplace equation. Hence in principle, if 
01 = is satisfied on all boundaries and provided that 02 vanishes at any one point, the two constraint equations 
are automatically satisfied in all places as long as the elliptic equations for T, R, and C are solved. 



C. Boundary conditions 



In order to solve the elliptic equations derived in the previous section, we must specify the boundary conditions. 
Let us first discuss boundary conditions on the symmetry axis of r = 0. The regularity at r = of Eqs. and 
requires 

R = C (at^ = 0), (13) 

which determines the boundary condition for R. The axial symmetry requires that the r derivative of the metric 
functions vanishes at r = 0. However, as long as we use the ^ coordinate, the finiteness of the derivative in this 
coordinate automatically guarantees the regularity because oc v^T^. Hence, we adopt a free boundary condition 
for T and C. The values of T and C at r = are evolved in the same way as the values at an ordinary grid point. 
The only difference is to use the one-sided differentiation to evaluate the first-order derivatives with respect to The 
second-order derivatives disappear at ^ = from the Laplace operator. Equation © at ^ = is trivial from Eq. H13|l . 
while Eq. (|10|l reduces to 



(at ^ = 0). (14) 



We can also show that this equation is guaranteed to be satisfied if R and C are solutions of Eqs. Q and Although 
it is not necessary in principle, we use this condition H14|) in addition to Eq. H13|l in order to improve the accuracy of 
the calculation. We do not solve the evolution equation for R at the ordinary grid points next to the axis in place of 
imposing this supplementary condition. 

On the horizon, the Killing vector dt must become null. Therefore, we have 

T-O (atp^ph), (15) 

In the present gauge, the horizon is given by a constant radius, and thus we have 

T^^T^^^Q {a.t p^ph). (16) 

Here we assume that metric functions and their derivatives are finite on the horizon. With the aid of Eqs. H15|l and 
()15|l . the surface gravity on the horizon is given as 

K^e-'^Tp (atp-Pft). (17) 

This condition can be used to impose a Dirichlet boundary condition for R on the horizon. For this purpose, we 
rewrite it as 

R = C\^=o + log ( jT^) (at P^Ph), (18) 
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using Eq. (O. 

Following the line of a standard proof can show that the zeroth law of black hole thermodynamics, d^n — 0, 

is valid, and it guarantees Eq. © to be satisfied on the horizon. Assuming the surface gravity is nonzero (nonextremal) , 
we find that Tp must be nonzero on the horizon; Tp ^ 0. Then, from the regularity of Eq. {Tj) or Eq. |(SJ) on the 



horizon, a condition for C p is derived as 



C,p = (atp = p,0- (19) 

On the other hand, the condition that the expansion is zero on the horizon gives [4^ 

■if 

{R + 2C),p + =0 (atp^ph), (20) 

PhZ 

which guarantees the constraint equation l|10() on the horizon. In the following numerical calculation, we basically 
impose these Neumann boundary conditions on the horizon for R and C. As a supplementary condition to improve 
the accuracy, we also use Eq. (fH^ . but we do not use Eq. (fT7|l . Hence, in our calculation the surface gravity is not 
given by hand as a parameter. We discuss this point in more detail at the end of this subsection. 

For the extremal case k — 0, the boundary conditions are not fully determined in the present way. From k ~ 0, 
i? = CX3 or Tp — are required on the horizon, and thus the boundary condition for C is not determined from the 
regularity on the horizon. It might be interesting to explore the extremal case separately, but we do not consider it 
in this paper. 

Israel's junction condition for the RS braneworld model gives K^^i, — —7^^/^, where 7^^ and K^^i, are the induced 
metric and extrinsic curvature on the brane, respectively. From this condition, we obtain the boundary conditions on 
the brane as 

^=d^R^d^C=-^{l-e'') (at e = (7r/2)2) . (21) 

The constraint equation @ is manifestly satisfied under these conditions. 

We must also specify asymptotic boundary conditions. For the asymptotic infinity, the boundary conditions to 
obtain asymptotically AdS spacetime are T 1 and R,C — s- 0. The metric functions must smoothly approach these 
asymptotic values. In actual numerical calculations, these asymptotic boundary conditions are imposed at a finite, 
but sufficiently far region, and we must check that the solutions are insensitive to the position. 

With these boundary conditions, we can solve the elliptic equations as a boundary value problem. Moreover, these 
boundary conditions guarantee the constraint equation Oi to be satisfied on all boundaries and 02 to be satisfied at 
least on the horizon. Hence as we explained, the two constraint equations are automatically satisfied as long as the 
elliptic equations are exactly solved. 

It should be noted that all elliptic equations, constraint equations, and boundary conditions can be rewritten in 
terms of nondimensional coordinates {p, x}, where p = without any dimensionful parameters. Then this system 
of equations, and hence a black hole solution, is characterized only by a single parameter 

L = e/ph. (22) 

In our numerical calculation, we take ph = 1 and change £ to specify this parameter, since we want to fix a coordinate 
region where numerical calculations are performed. The variation of £ keeping L fixed corresponds to a rescale of the 
length scale as we notice from the metric which is rewritten as 

ds^ =£^dsl{ip,x), (23) 

where ds\ is the nondimensional part of the line element written by {t = t/£,p,x}, and its metric function is given 
by a black hole solution T, R, and C with the parameter L. After specifying a black hole solution by this method, one 
may want to transform the solutions to those in the coordinates where L is specified by taking £ = 1 and changing 
Ph- For this purpose, it is enough to multiply £~^ to a length scale in the coordinates where p/j = 1 is taken. 

A question may arise. Although we say that we use L = £/ ph to specify a solution, ph does not have a clear 
physical meaning. One may think that it is more appropriate to specify the value of the surface gravity k, instead 
of ph- Here we should recall the argument given at the end of the preceding subsection that p/ pt and hence p/ Ph 
cannot be changed arbitrarily by using the residual gauge degrees of freedom. If we change the value of ph by using 
this residual gauge degree of freedom, the value of p at infinity also scales correspondingly. Then, we will see that the 
asymptotic boundary conditions T — > 1 and i?, C — > are not satisfied after the gauge transformation. Hence, for a 
given solution we do not have a degree of freedom to change the value of ph arbitrarily. Thus, if we specify a solution 
by fixing both k, and we are not allowed to set p/j = 1 any further. On the contrary, as we wish to set p/j = 1 in the 
actual numerical computation, we should not use the condition H17|) that specifies the value of k directly. 
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III. NUMERICAL EXAMPLES OF SMALL BLACK HOLES 
A. Numerical examples 

The formulation outlined in the previous section indeed works well to solve the Einstein equations as a boundary 
value problem, using an iterative convergence scheme of the relaxation method. A simple check of an algorithm is to 
calculate a very large limit of L for fixed p/i, and compare the result with the 5D Schwarzschild black hole in isotropic 
coordinates which is a solution in the small horizon limit L oo — > cx)) taking z/i 1 



Tip) = 



P' - Pi 
P^ + pI ' 



R{p) = C{p) = log + ^) • (24) 

In this limit, the tension of the brane vanishes and the background spacctimc is flat. We have used this metric as 
an initial guess configuration with which we start the relaxation. Namely, we first calculated a solution for relatively 
large L, and then proceeded to smaller L by using the result relaxed for the previous value of L as the new initial 
configuration. Convergences of numerical calculation become worse as L becomes smaller, and errors of constraint 
equations also grow if we keep the same resolution. It becomes in general more difficult to find solutions as the 
nonlinear effects of differential equations grow, and we could not keep the convergence of the calculation for L < 1. 
Thus in this paper we show only small black holes {L > 1), although it is important to improve our scheme by 
identifying and removing all numerical instabilities. The number of grid points {^i , pj ) used in our calculation is 
100 X 1000. The grid point pj is taken to be a geometric progression pj = ph + e(l — 7-')/(l — 7). The parameters e 
and 7 are determined by requiring that the ratio of the last grid resolution to the first one is 6pina.x/^Pi ~ 2.71 and 
the asymptotic boundary is set at Pmax = 85. It is crucial to test the sensitivity of the solutions to the physical size 
of the lattice. This will be discussed later (see Fig. IHJ. 

To orient the reader, we show the result of T, R, and C for L = 15 in Fig. as a typical example. The center of 
the black hole is at r = and z — £ = 0. The interior of the horizon p < ph is outside the region of our numerical 
calculation. The contours of T are almost spherical for this value of L as in the case of the limit L ^ 00 in spite of 
the nontrivial boundary conditions on the brane, while the contours of R and C manifestly deviate from spherical 
shape. 

In Fig. 121 the solution for i = 30 is displayed in the polar coordinates that are more appropriate to see the angular 
dependence of the metric functions. Figure displays the results for L = 10 in the polar coordinates. The deviations 
from spherical shape for L = 10 are enhanced more than those for L = 30. We have performed numerical calculations 
for the parameter region L = 3 ~ 500. The results for some L are listed in Tabled 

To consider the shape of the small localized black hole in the bulk, we calculated geometrical quantities for the 
numerically obtained black hole solutions. The area of the black hole horizon is given by 

(-/^)^ sin' Vi ( f\\n+2C 



2V? V 



A5 = 2pi / dn / df-^^-^ - e^^+^^. (25) 



Here the factor 2 is due to Z2 symmetry. The proper area of the intersection between the horizon and the brane is 
given by 



A4=ATTpie"- (x = V2). (26) 
Figure 0] shows the ratio of a mean radius in four-dimension y/A^ (on the brane) to that in five-dimension Ag^"^ 



^5 

The figure indicates that the black hole tends to flatten as its horizon radius increases. This geometrical behavior is 
expected from the AdS geometry of whole spacetime. The so-called warp factor £z~-^ of the AdS geometry ^ will in 
general cause the localized black hole to be flattened, and this effect becomes significant when the horizon radius is 
large. In fact, the EHM solution exhibits a similar behavior to Fig. 0] 

B. Constraint equations and checks of calculation 



In this section we discuss the validity of our numerical calculations. We have calculated the constraint equations 
for each numerical solution. Averaged absolute errors of the elliptic equations and the constraint equations are listed 



7 



in Table ni In the table, we have introduced a norm N{^,p) of the constraint equations since we need appropriate 
references to compare with the violations of constraint equations. The norm is defined by 



where n^'^^ represents the terms in <di or Q2 at {£,,p)- The comparison between this norm and the violation of the 
constraint equation gives a degree of cancellation between different terms. We see in Table n that averaged constraint 
violations against averaged norms are small, and thus the constraint equations appear to be well satisfied. Moreover 
by using this norm, we can obtain relative accuracies of the constraint equations, |6i|/A^i and |62|/-^2- Figure El 
shows that the absolute errors of the constraint equations, |9i| and I62I, are observed mainly around the origin of 
polar coordinates and near the horizon where the elliptic equations have terms whose coefficients behave like —^ 00 
or 1/T — > 00. On the other hand, the figure represents that the relative accuracy around the horizon is not so bad, 
although the relative accuracy becomes worse near the axis. This violation of the constraint equations is a general 
feature of the present calculations and it becomes larger as L becomes smaller. The lack of accuracy near the axis is 
a common problem in axisymmetric problems. 

We have also checked second-order convergences of the thermodynamic quantities for L — 20 keeping pmax — 85. 
As a function of S^, we confirmed the quadratic fits for k and A5. The similar quadratic fit is possible for the variation 
of p resolution, in which we need to change the lattice resolution keeping the ratio Sp^i^-^/Spi fixed (« 2.71). The 
numbers of grid points {£,i,Pj) used in these checks are (80 — 200) x 1000 for SS, variation, and 100 x (750 — 2000) 
for 6p variation, which correspond to 6^ = 0.01 — 0.03 and dp = 0.02 — 0.07, respectively. Note that the absolute 
errors decrease as the lattice resolutions increase. However, clear second-order convergence has not been observed 
because we have not used equal grid spacing for pj. If we use equal grid spacing, second-order convergence is actually 
observed. We have also checked the insensitivity of our numerical solutions on the finite position of the asymptotic 
boundary. To confirm it, we vary the position pmax for L = 10 and L = 20, and check the robustness of thermodynamic 
quantities. For each variation of pmax, we have changed the grid resolution for p keeping the ratio Spmax/Spi ~ 2.71 
with Spi ~ 0.05. The thermodynamic quantities are stable for each value of Pmax- Figure shows a result for a 
nondimensional combination. One sees that the combination of thermodynamic quantities is stable and the variation 
is no more than 1% for the parameters. 



Black hole solutions that represent a black hole localized on the brane in the Randall- Sundrum infinite braneworld 
model have not been found with appropriate boundary conditions in the literature. We have explored this problem in 
this paper. We have performed the numerical calculations and found nontrivial localized black hole solutions whose 
horizon radii ph are small compared to the bulk curvature scale £. More precisely, the small localized black hole 
solutions are approximately constructed by numerical calculations under appropriate boundary conditions. From the 
numerically obtained solutions, we could observe the transition of the shape of the small localized black hole in the 
bulk (Fig. 2)). The black hole tends to flatten as its horizon radius increases. This is expected due to the nature 
of the AdS geometry, and is also observed in the EHM solution of a localized black hole in the 4D braneworld 39l| . 
Although the EHM solution is the solution in lower dimensions, it is interesting to compare our solution with it in 
more details. 

Although the method we have developed works well for the small localized black holes, we could not succeed in 
finding black hole solutions with large horizon radius. The lack of convergence in the relaxation scheme is due to the 
appearance of unstable modes. For example, an unstable mode that appears near the axis is a common problem in 
axisymmetric code (see, e.g., 0). In particular, the first terms in Eqs. Q) and JHJ) that have the factor 1/r^ give severe 
contributions near the axis, and they invoke numerical instability. The same problem appeared in the calculation of 
a relativistic star in the braneworld model ,44| . Thus improvement of the numerical scheme and formulation to avoid 
this instability is an important issue. It may be worth pointing out that since we do not have any uniqueness theorem 
in the context of the braneworld |53| , there is a possibility of finding a sequence of solutions other than that found 
in this paper if we start with a completely different initial guess. 

It is worthwhile to point out that the boundary conditions that we have imposed do not a priori guarantee the 
absence of naked singularities in the bulk or on the horizon. If any solution obtained in our formulation cannot avoid 
a naked singularity in the bulk or on the horizon, it means that there is no physically acceptable static black hole in 
the RS infinite braneworld. Such a possibility was pointed out for large localized black holes 

As a direction of future work, it is crucially important to find large black hole solutions if they exist 0iE3|i and 
it needs further developments of technique and investigation. It is also interesting to apply this method to find the 




(27) 



k 



IV. 



SUMMARY 
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black holes localized on the TeV brane [5j. The stabilities of a higher-dimensional black hole including our solution 
are also important. We hope the method used in this paper will become the basis for further studies on black holes 
in higher dimensions. 
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"^ 3.7x10" 


0.37 129 53 
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(1.2x10" 


"^ 2.7x10" 
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L = oo 
















0.50 158 50 



TABLE I; Table summarizing the data of numerical calculations for several values of L = l/ph- In the calculations we have 
varied £ keeping the horizon radius ph = 1- Thus it should be noted that the dimensionful parameters, i.e., surface gravity k, 
5D area A5, and 4D area A4,, must be rescaled for comparison, as we have explained below Eq. 12211 . In the table, L = 00 
represents the 5D Schwarzschild black hole, which is a solution in this limit. {IO2I), and {\TRC\} are averaged violations 

of the constraint equations and the elliptic equations. {\TRC\) is the mean absolute value for the equations of T, R, and C with 
equal weight for each grid point. In taking these averages, the first 5% of lattice points in ^ and the region near the horizon 
p < 2ph are excluded since they receive unphysical enhancements of errors due to singular terms. We introduce (A'^i) and 
as the averaged value of TV [Eq. 1271 1 for Bi and Q2, respectively, in the restricted grid region mentioned above. The number 
of grid points {^i,Pj) used in the calculations presented in this paper is 100 x 1000. The coordinate position of grid point pj 
is taken to be a geometric progression pj = pjj + e(l — 7-')/(l — 7). The grid spaces used in the calculation are = 0.025, 
= 0.049, 5pmax/5pi = 2.714. The position of the asymptotic boundary is set at pmax = 85. 
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FIG. 1: An illustration of the metric functions T, R, and C for L = 15 in the coordinates {r,z}. This is a typical example 
of our numerical solutions for small black holes {L > 5). In these figures the center of the black hole is placed at the bottom 
left corner (r = and z — £ = 0), and the brane is at 2: — £ = 0. The numerical calculation has been performed in the polar 
coordinates {p,^}, and hence the inner region of the black hole p < ph is outside the area of our computation. 
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FIG. 4; An illustration of y A^jA^ . This plot shows a degree of deformation of a black hole. One sees that as L decreases 
the ratio of a mean radius in four-dimension A^^^"^ (on the brane) to that in five-dimension A^'^j^ increases. It indicates that the 
black hole tends to flatten as its horizon radius increases. The dashed line shows the same quantity for the 5D Schwarzschild 
black hole. Note that L is expressed in natural logarithm. 
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FIG. 5: An illustration of the constraint equations for L = 15. The figures that plot the absolute values of two constraint 
equations, \ and 1 02 1 , show the absolute errors of the constraint equations. The other two figures plotting 1 0i (p, x) | /Ni (p, x) 
and |©2(p, x)\/N2(p, x) show the relative accuracy of the constraint equations since A''i and A'2 are norms of respective constraint 
equations (see Table EJ. As expected, the absolute errors of the constraint equations are observed mainly near the axis of polar 
coordinates and around the horizon. However, the relative accuracy is not significantly low around the horizon, but it is worse 
near the axis and in the asymptotic region. 
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FIG. 6: A plot showing the sensitivity of the numerical solutions to the position of asymptotic boundary. The calculations 
are performed for L = 10 and L — 20. The position of asymptotic boundary pmax is changed from pmax = 50 to pmax = 115 
keeping the ratio of the resolution fixed at 5pmax/5pi ~ 2.71 with Spi ~ 0.05. We see that the variation of the thermodynamic 
quantity is no more than 1%, and is stable for this parameter range of pmax. 



